Mon. Not. R. Astron. Soc. 000, 1-12 (2007) Printed 25 March 2008 (MN I^TeX style file vl.4) 



Two-integral distribution functions for axisymmetric 
stellar systems with separable densities 

Zhenglu Jiang 1 * and Leonid Ossipkov 2 t 

1 Department of Mathematics, Zhongshan University, Guangzhou 510275, China 

2 Saint Petersburg State University, Staryj Peterhof, Saint Petersburg 198504, Russia 

Accepted 2007 September 25. Received 2007 September 15; in original form 2007 July 14 

ABSTRACT 

We show different expressions of distribution functions (DFs) which depend only on 
the two classical integrals of the energy and the magnitude of the angular momentum 
with respect to the axis of symmetry for stellar systems with known axisymmetric 
densities. The density of the system is required to be a product of functions separable 
in the potential and the radial coordinate, where the functions of the radial coordinate 
are powers of a sum of a square of the radial coordinate and its unit scale. The even 
part of the two- integral DF corresponding to this type of density is in turn a sum or an 
infinite series of products of functions of the energy and of the magnitude of the angular 
momentum about the axis of symmetry. A similar expression of its odd part can be 
also obtained under the assumption of the rotation laws. It can be further shown that 
these expressions are in fact equivalent to those obtained by using Hunter and Qian's 
contour integral formulae for the system. This method is generally computationally 
preferable to the contour integral method. Two examples are given to obtain the even 
and odd parts of their two-integral DFs. One is for the prolate Jaffe model and the 
other for the prolate Plummer model. 

It can be also found that the Hunter-Qian contour integral formulae of the two- 
integral even DF for axisymmetric systems can be recovered by use of the Laplace- 
Mellin integral transformation originally developed by Dejonghe. 
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1 INTRODUCTION 

This paper is a continuity of our recent work (Jiang & Ossipkov 2007a,b) to construct a self-consistent stellar system by 
means of finding a two-integral distribution function (hereafter DF) for a stellar system with a known gravitational potential. 
Jiang and Ossipkov (2007a) have shown a method of finding anisotropic DFs for spherical galaxies. This is an combination of 
Eddington's (1916) formula and Fricke's (1952) expansion idea. Of course, they can be also regarded as simply an extension 
of the idea of Eddington. Some similar formulae are given by Jiang and Ossipkov (2007b) for finding two-integral DFs for 
axisymmetric stellar systems. 

There is a long history of finding the DFs for a stellar system with a known gravitational potential. Eddington's (1916) 
formula is for spherical galaxies and is in fact a solution to the Abel integral equation (see Binney & Tremaine 1987, hereafter 
BT; Dejonghe 1986). Fricke (1952) gave an expansion idea, that is, DFs which are products of the two powers of the energy 
and the square of the angular momentum about the axis of symmetry correspond to densities which are proportional to 
products of the potential and the radial coordinate for axisymmetric systems. After that, different integral transformation 
techniques are used to obtain the DFs of spherical or axisymmetric stellar systems (e.g. Lynden-Bell 1962a; Hunter 1975; 
Kalnajs 1976; Dejonghe 1986). Hunter and Qian (1993) gave complex contour integral formulae to get the two-integral DFs 
for axisymmetric systems. 
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Therefore, in this paper, we present different expressions of the two- integral DF in the stellar system with the axisymmetric 
density that is required to be expressed as a product of two functions: one is a function of the potential and the other a power 
of a sum of a square of the radial coordinate and its unit scale. We first introduce the fundamental integral equations of 
the problem of finding the two-integral DFs in Section 2 and then present our method of solution in Section 3. We show in 
Section 3.1 that the Hunter-Qian complex contour integral formulae can be recovered by use of the Laplace-Mellin integral 
transformation defined by Dejonghe (1986) and that the two- integral even DFs for stellar systems with the known axisymmetric 
densities of this type mentioned above can be expressed as a finite sum or an infinite series of products of functions of the 
energy and of the magnitude of the angular momentum about the axis of symmetry. In Section 3.2 we give the similar 
formulae of the two-integral odd DFs under the assumption of the rotation laws. In Section 4 we use this device to obtain 
the two-integral DFs for both the prolate Jaffe and Plummer models. These DFs are indeed in accordance with those derived 
from the Hunter-Qian contour integral. Section 5 summarizes our results and exhibits our conclusions. 



2 THE FUNDAMENTAL INTEGRAL EQUATIONS 



Assume that a stellar system has a gravitational potential with a upper bound. Let $ and E be, respectively, the potential 
and the energy of a star in the system. As in BT, we can choose a constant $o such that the system has only stars of the 
energy E < <f>o, and then use the usually so-called relative potential ip — — <E> + $ and energy e = — E + $o to describe the 
stellar system. Here, e — is a relative energy of escape from the system. Given a relative potential ip — tp(r) of a stellar 
system its mass density p = p(r) can be derived from Poisson's equation. Then the problem of finding the DF / = /(r, v) of 
the system is that of solving the following integral equation 



/ 



(1) 



where r is a position vector, v is a velocity vector. In an axisymmetric system, we generally use the cylindrical polar coordinates 
(R, <p, z) with the z-axis being that of symmetry. Let v = (vr, v^,,Vz) and L z be, respectively, the velocity and the component 
of angular momentum about the z-axis. Then we can know that the relative energy e and the z-axis angular momentum 
L z = Rv<f, are two isolating integrals for any orbit in the axisymmetric system. Thus, by the Jeans theorem, the DF of a 
steady-state stellar system in an axisymmetric potential can be regarded as a non- negative function of e and L z , denoted by 
f(e,L z ), and then for an axisymmetric system, (1) can be rewritten as 



/ 



f{e,L z )d\. 



(2) 



Put v„ 



— \J v\ + vl and define cylindrical coordinates (v m ,v<f,,9) in velocity space by the following relations: vr — v m cos 9, 
v z = v m sin 9. It thus follows that d 3 v = v m dv m dv<j,d9. Since e = ip — (l/2)(v R + v| + v z ), the integrand of (2) is independent 
of 9 so that we integrate out 9. We then change the variables (v m ,v^) to (e,L z ). It is easy to see that the integral in dL z is 

e) to R^U^p 



actually an integral of / from —R^j2(ip — e) to R^2(ip — e). We finally restrict ourselves to only positive values of L 
changing the lower integration limit to and integrating the function 2/+(e, L z ), thus getting 



P = 



4tv 
~R 



R^2(^-s) 



f+(e,L z )dL z 



de 



by 



(3) 



for the system having only stars with e > [that is, f(e,L z ) = for e < 0], where f+(e,L z ) — [f(e,L z ) + f(e,—L z )]/2. 
This implies that a given density determines /+(e, L z ) which is just the part of the DF that is even in L z . Hence f+(e, L z ) is 
usually called the even DF. 

Using the similar derivation to that of (3), we can rechange p{v<p) — J w^/(e, L z )d?\ as 



4-7T 



Lzf.(e, L z )dL z 



de, 



(4) 



where f_(e, L z ) = [f(e, L z ) — f(e, —L z )]/2. Usually, f_(e, L z ) is called the odd DF of the stellar systems. Equation (4) means 
that the odd DF can be derived from the rotational velocity. We can also find that if p and f+(e, L z ) in (3) are replaced by 
pR(v^) and L z f_(e, L z ), respectively, then (3) becomes (4) for L z f_(e, L z ). 

Similarly, put Q = e — L 2 Z /(2R^), where R a is a scaling radius, and assume that the DF is dependent on Q and L z , 
denoted by f(Q,L z ), and that the system has only stars with Q > 0, or equivalently, / = for Q < 0. Obviously, Q — > e as 
R a — ► oo. Then, for an axisymmetric system, the integral equation (1) can be rechanged as 



-/ 



f(Q,L z )d\. 



(5) 
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By changing the variables of the integral in (5), it follows that 

dQ, 



4% 
~R 



f+(Q,L z )dL z 



where f+(Q, L z ) = [f(Q, L z ) + f(Q, -L z )]/2. Naturally, f+(Q, L z 
[f(Q, L z ) - f(Q, -L z )]/2 of the DF satisfies 

r fl A / 2 M'-Q)/(i+« 2 /«J) 

L z f_(Q,L z )dL z 



4tt r r 

R 2 Jo Jo 



dQ. 



(6) 



is the even part of f(Q, L z ). Also, the odd part f_(Q, L z ) — 



(7) 



Obviously, (5), (6) and (7) coincide with (2), (3) and (4), respectively. In other words, (2), (3) and (4) are, respectively, 
limits of (5), (6) and (7) when R a — » oo. Thus it suffices to consider the case of f(Q, L z ) in the next section. 

Once f+(Q,L z ) is known, f(Q,L z ) can be obtained under some assumptions (e.g. Lynden-Bell 1962a; Miyamoto 1971, 
1974; Ossipkov 1978; Dejonghe 1986; Kutuzov 1995) and further p can be easily calculated by integration and tp by solving 
Poisson's equation for the axisymmetric system. 

The inverse problem that is now investigated for any axisymmetric system is how to derive the two-integral even DF 
f+(Q,L z ) from the density p and how to deduce the two-integral odd DF f_(Q,L z ) under the assumption of the rotation 
laws. 



3 THE TWO-INTEGRAL DFS 

In this section, we shall show different expressions of the even and odd parts of the two-integral DFs corresponding to three 
types of the mass densities given by 

p(i>,R) =p(*p)/(l + R 2 /Rl), (8) 
p(^,R)=pM/(l + R 2 /R 2 a ) 1/2 , (9) 
p(^,R)=pW/(l + R 2 /R 2 a ) 3/2 . (10) 
We shall also give the similar expressions of the even and odd DFs for the spherical mass density p = p(tp). 

3.1 The even DFs 

We can know that under some suitable assumptions a complex integral expression of the two-integral DFs corresponding to 
the above three types of the axisymmetric mass densities can be obtained by use of Hunter and Qian's (1993) complex contour 
formula. It can be further known that this contour integral is the same as that derived from Dejonghe's (1986) Laplace-Mellin 
integral transformation for the axisymmetric mass density of the separable form p(ip,R) = F(ip)g(R), where F(ip) and g(R) 
are two functions such that both the Laplace and Mellin integral transformations of F(ip) and g(R) are valid, respectively. 
To show this, we first have to recall the following definitions of the Laplace-Mellin integral transformations 



(11) 



-|-oo 

ai, dip J R p - 1 p(^,R)dR (12) 



and 

£,^oWr^p{p} = J e 



with their inversions 

ao+ioo /3 + ioo 

/ + ( £ ,L,) = -JL J e «z da J L-P£^ a M L ^ {f + }dp (13) 

cto — ioc Pq— ioo 

and 

QO+ioo /3 + ioo 

p(ij,R) = - 1 L J e ^ da J R-^^jm R ^{p}dP, (14) 



© 2007 RAS, MNRAS 000, 1-12 



4 Z. Jiang and L. Ossipkov 

A Im(</>) 



C(e) 




Re(i/>) 



Figure 1. The contour C(e) of the complex integral given by (19). Both the branch point and the pole of its integrand are the relative 
energy e in the real V-axis to the left side of the maximal relative potential No- 



where, Qo and /3o are two suitable real numbers, £ and 9JI denote the Laplace and Mellin integral transformations, respectively. 
Here and below, i is an imaginary unit. Obviously, the Laplace- Mellin integral transformations of (8), (9) and (10) are valid 
for Re(/3) < 1 when p(tp), p(ip) and p(ip) are suitable for the Laplace transformation. 

We can then know that the Laplace-Mellin integral transformations of the mass density p = p(tp,R) and its DF /+ = 
f+(e,L z ) have the following relation [see Dejonghe 1986, Page 273, (2.3.3)] 



a {Z-0)/2 2 0/2 



(15) 



where T(x) denotes the Gamma function, — 7r/2 < arg(a) < 7t/2 and Re(/3) < 1. Note that the Laplace-Mellin integral 
transformation of the separable density p(ip,R) = F(tp)g(R) is given by £^,_ >a £!Hj{_ > /3{p} = £ii,^ a {F}9JlR^f3{g}. Thus the 
Mellin integral transformation of /+ = f+(s,L z ) can be obtained as 

3Kw{/f} = ^aXrff)^ -- {« (3_/3)/2£ ^«{^} • ( 16 ) 

where £q^ £ denotes the inversion of the Laplace integral transformation. It can be shown that the Laplace inversion trans- 
formation on the right side of (16) has the following form (see Dejonghe 1986, Page 306, Appendix Al, Theorem 1) 

1 



£-i £ {a< 3 -«/%_{F}} = 



r(±±£) 



de 2 



F(^){e-iP) (p - 1)/2 d^, 



where 1 > Re(/3) > -1. Assume that F(0) = F'(0) = 0. Then (17) gives 



£-i £ {«< a -«/ 2 £,_{F}} 



r(i±^)7 W 



d 2 F 



x^W(e-V<) 



(/9-l)/2 



dip. 



(17) 



(18) 



Assume without loss of generality that F(ip) has an analytic continuation near the real ?/>-axis in the complex ?/>-plane. Then 
equation (18) can be rewritten as 



£^ e K^ )/2 £^ Q {F}} = - 



2cos(^)r(i±£) 



f £1. 

Jc(e) 



W(e-V-) 



(/3-l)/2 



dtp, 



(19) 



where and below, C(e) is a suitable contour of the complex integral (19) with respect to the complex variable tp; the contour 
is a loop which starts from the lower side of the real ip-axis at ip = 0, passing through a physically considered window V in 
the real ?/>-axis to the right of tp = e, to the upper side of the real i/>-axis at i/j = 0; the window V physically considered must 
lie between the relative energy e and the maximal relative potential ipo. Figure 1 shows one cut (broken line), one window V 
and one contour C(e). Inserting (19) into (16), we can get 



3rti.^8{/+} = - 



4\/27r 



C(e) 



d 2 F 

dip 2 



V^? r(!=«)r(^)™(fe) 



(20) 



'1-/31 



By using the identity equality that T 

Mellin inverse integral transformation of equation (20) can be expressed as 



cos (^) = 7T (see Gradshteyn & Ryzhik 1965, Page 937, 8.334, 2), the 



(21) 
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where yjl fj ^ L denotes the inversion of the Mellin integral transformation. (21) can be equivalently rewritten as 
f+(e,L z ) = ^- [ ^r(i>)- F ^=M- 1 r- {m R ^{g}e- iB ^\ dip 



'C(e) 

(b)_ 

4^271 



fcts'^'t/s)^ (22) 



where (a) results from the definition of the Mellin integral transformation and (b) is obtained by virtual of the fact that 
OT r ^p {g(re w )} = e' ipe m r ^p{g(r)} [see Davies 2002, Page 203, (12.12)]. Hence (22) is in accordance with the Hunter-Qian 
complex contour integral (Hunter & Qian 1993) when the mass density p(if),R) = F(ip)g(R) is any one of the three densities 
defined by (8), (9) and (10). It can be similarly shown that the two- integral even DF 

'♦<■■« --i^/ m &(W*Fw);£7 

corresponds to the general axisymmetric mass density p = p(tp, R) having an analytic continuation in the complex plane with 
respect to the variable ip and satisfying that p(0,R) = [dp{ip, R)/dtp]^=o = 0. Obviously, (23) is the same as the complex 
contour integral formulae mentioned above. 

Let p and v be two constants. Assume that — p — v is not a natural number and that p > 3/2. If the two functions F and 
g in (16) are defined as F(ip) = tp p and g(R) = (R 2 /R 2 a )^/(1 + R 2 / R 2 a Y +v , respectively, then £^ a {F} = r(p+l)/a p+1 and 
Wl R ^i3{g} = {R 2 a f /2 r(p + (3/2)T(v - P/2)/[2V{p + v)]. Thus (16) becomes 

m If , r(p + i) r( / x + |)r(^-|)(2ig)^ 2 (1 _, )/2 _ n 

Y{ P + l)e^' 2 v(p + ^v(v^^){2eR 2 a f' 2 

TVHTfr + v) r(i^)r( P -i + |) • 1 ' 

Then (24) can be equivalently rewritten as 

r( P + i) £ ^ i °T r( M +g)r(„-f) /^V"' 

/+(£ '^ j 2^nT(p + u)2ni J r (±^) T (p - ± + f ) ^7 * P 

/3 -ioo 



2 3 / 2 7rr( M + ^) ^' " 2' 2' 2ei?i 
where H(/i, z/, c, d; x) is defined by 

with a contour C determined by /3o such that — p are on its left side and v on its right side. In the case that p + r] and v + £ 
are not negative integers, the complex integral (26) can be calculated and expressed as follows. When < x < 1, if p — £ is a 
nonnegative integer, then H(/i, v,{;,r);x) — 0, or else H(/_t, z/, £, 77; cc) = x M 2-Fi(/i + f, — £; M + £)r(/^ + f )/[r(£ — p)T(p + rj)]\ 
when a; > 1, if v — 77 is a nonnegative integer, then M(p, v, £, 77; x) = 0, or else H(/_t, 1/, £, 77; x) = x~ v iF\{p + i>, 1 + 1/ — r\\v + 
i\ l/a;)r(^ + i/)/|r(»7-i/)r(i/ + |)]. Here, 2 F\ is a hypergeometric function. In particular, B(0, l,p-l/2, 1/2; x) = 2F1 (1,3/2 - 
p; 1/2; x)/[r(p- l/2)r(l/2)] for any x € [0, 1), and H(0, l,p - 1/2, 1/2; a) = aT^iM 1 , 3/2; p+ 1/2; l/z)/[r(p + l/2)r(-l/2)] 
for any x € (1, +00). 

We shall below give different expressions of the two-integral DFs for the systems. First, let us consider the two-integral 
even DFs corresponding to the mass density defined by (8). Assume that p(ip) in (8) can be expressed as a convergent power 
series of the form p(ip) = &ntp n for p > 3/2. Then the mass density defined by (8) can be rewritten as 

P ^,R) = j2^ n /( 1+R2 / R ^- ( 2? ) 

n—p 

By (25), it can be found that the even DF corresponding to the mass density (27) is 

, ^ a n r(n + l)e^ 2 ^ ( 1 1 L 2 

M £ > L ^ = L 237% [°> n 2' 2 ; 27R 2 
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= < 



n— p 



-i . , (28) 

2(27r)3/ 2 r(n+l/2) V 2lflJ < 2 " 1 V 2' " 1 2> ~TJ 



It is worth mentioning that the series in (28) are convergent to the corresponding contour integrals obtained by virtual 
of the Hunter-Qian contour integral formulae for the systems since all the terms in this series can be equivalently rewritten 
as the corresponding Hunter-Qian contour integrals. If the convergence of the above power series is fast enough, then one 
can adopt a sum of finite previous terms of the series (28) as a viable approximation of the even DF. Note that it is easy to 
calculate numerically these hypergeometric functions in (28). Hence this method is in general computationally preferable to 
the contour integral method 

By using the method given by Jiang and Ossipkov (2007b), it can be easily known that the stellar system with the mass 
density of the form (9) has the two-integral even DF f+(Q) of the form 



since p(ip) vanishes at ip = for the stellar system. 

Then we shall show the two-integral even DF corresponding to the mass density (10). Note that (10) can be rechanged as 

P (i>,R) = pW[Rl/(i + R 2 /Rl) 1/2 - R 2 /a + R 2 /Rlf /2 ]/Rl (30) 

Assume that [<P p(ip) /dip J ]^=o = for j = 0, 1. Similarly, it then follows that the two-integral even DF 
j2 rQ ;t/„,,\j„,. r2 j3 rQ 



f (n T ) = — * f W PW# Ll d 3 T 

corresponds to the mass density (30) (or say, (10)). 

It can be also proven that (29) and (31) are in accordance with those obtained by use of Hunter and Qian's (1993) contour 
integral formula and that they can be derived using Dejonghe's (1986) Laplace-Mellin integral transformation. Obviously, these 
real integrals are easier to calculate than the corresponding complex contour integrals. 

We can finally mention that for the spherical stellar system with the mass density p = p(tp), there is a unique isotropic 
DF of the Eddington (1916) type 

i d 2 r pwdv> 

f+{£) -7^d^J 7^ (3) 

and an infinity of different anisotropic DFs of the Ossipkov-Merritt type corresponding to an infinity of various values of an 
anisotropy radius (Ossipkov 1979; Merritt 1985). 



3.2 The odd DFs 

Note that if p and f+(Q, L z ) in (6) are replaced by pR{v<p) and L z f_(Q, L z ), respectively, then (6) becomes (7) for L z f_(Q, L z ). 
Then we can use Hunter and Qian's contour integral (or Dejonghe's Laplace-Mellin transformation, mentioned above) to get 
the odd part of the two-integral DFs corresponding to the three types of the axisymmetric mass densities. We shall below 
show other expressions of the two-integral odd DFs. 

Assume that the rotational velocity of the stellar system obeys the rotation law: (v^) = R 2k /(l + R 2 /R 2 ) k , where k 
is a positive constant. In the case of k being a natural number, this law has been used to get the odd DFs of Binney's 
(BT) logarithmic model (Jiang & Ossipkov 2007b). When k = 1, it is one of the rotation laws considered by Evans (1993). 
Define [k] = k — ko, where [k] represents a nonnegative integer and ko G [0, 1). Suppose that p(tl>) in (8) satisfies that 
[d j p(tp)/dtp j ]^=o = for j = 0, 1, • • • , [k] + 1 and that the system has only stars with Q > 0. Take f_(Q, L z ) = sgn(L z )L 2 z k h(Q) 
as a two-integral odd DF for the mass density defined by (8). Then inserting (8) into the integral equation (7) gives 

R2 "M) - ^ R2k C h(Q) ^ _ Q) k+1 dQ (33) 

{l+w/Rzy+i - (fc+i)(i+i?v^) fc+i y 

It can thus be found from (33) that 
(k + l)J 

Using the Laplace integral transformation (Dejonghe 1986), we can get the following inverse formula of (34): 

h(Q ) = k + 1 ^ [ Q (35) 

lVj 47r2 fc r(fc + 2)r(l - ko) dQW+ 3 J (Q-V) fe "' 



/ 5 W = 77 Z ?Tt/ M<2)W> - Q) k+1 dQ. (34) 
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Therefore we can know from (35) that under the above assumptions the two-integral odd DF is given by 
f(n r)- (k+l>gn(L z )Ll k rfM+ 3 f Q pW)dnj> 

/-w, *") 47r2 fer(fc + 2)r(i - fco) rfQW+ 3 J (Q - V) fe ° { ' 

for the mass density (8). In particular, when k is a natural number, the two-integral odd DF (36) for the mass density (8) 
can be rewritten as 

sgn(L z )L? d k+2 p(Q) 

on the condition that [dp p(tp) / '<ft/> J ']^,=o = for j = 0, 1, • • • , k + 2. 

Next, we consider another expression of the two-integral odd DF corresponding to the mass density of the form 

+ OC 

p{i>,R) = ^a n ^ n /{l + R 2 /Rly<. (38) 

n—p 

for any given nonnegative constant 7. Note that if p and f+(e, L z ) in (3) are replaced by pR{v^,) and L z f_(e, L z ), respectively, 
then (3) becomes (4) for L z f_(e, L z ). Hence, by (25), we can find that the two-integral odd DF for the mass density (38) is 

.... ^ a n R™ +1 T{n + l)e"-W* m {L z )\L z \-i ( 1 1 11. L 2 Z \ 

n—p 

+ °° / r 2 \ 

^ a „n! E "- 2 H(fc, 7 ,n-l,l;^ (39) 



47rr(fc + 7) 

under the above assumption of the rotation law. The above last equality is obtained by using the following identity equality 
£ - 3 / 2 |L z |- 1 e (mZ^JZ) =e"- 2 (2R 2 a )- 1 / 2 M !,„+!,£_ - „- . 



2eR 2 J v a ' \y 2' 2'^ 2' ' 2 ' 2eRl 

In (39), the functions H can be expressed in terms of hypergeometric functions and so they are easy to calculate numerically. 
This expression (39) is suitable for the mass densities defined by (8)- (10), even for the spherical mass density p = p(ip)- 

It can be finally found that all the above series expressions of the odd DF are also convergent to the corresponding 
complex contour integrals obtained by using Hunter and Qian's (1993) contour integral formulae. If the convergence of the 
series (38) is fast enough, then we can use a sum of finite previous terms of (39) as a good approximation of the odd DF. 



4 APPLICATION TO THE PROLATE MODELS 

In this section, we shall apply the above series expressions into getting the even and odd parts of the DFs for both the prolate 
Jaffe and Plummer models. 



4.1 The prolate Jaffe model 

The gravitational potential of the prolate Jaffe model is as follows: 

*(*, K ) = ^lJ + + + ) (40) 

r J V VWR 2 + a ' 2 + b ) 2 + z2 + r J J 

Here, and everywhere below, a, b, M and rj are positive constants and G is the gravitational constant. By using Poisson's 
equation V 2< £ = 4nGp, it is easy to show from (40) that 

, D . M r 2 T 2 (X [i + bX 2 + a 2 b) + brjr [{ (X 2 +a 2 ) + (3r + 2rj)a 2 rjX(X + b) 2 

p{R > z) = ^-j ^{T + rjyx* ' (41) 



where and below, r = \/~{X + b) 2 + z 2 and X = y/R 2 + a 2 . This density-potential pair was given by Jiang and Moss (2002) 
and a two-integral even DF corresponding to this pair was shown by calculating the Hunter-Qian contour integral. We can 
below express the two-integral even DF as a series of hypergeometric functions, and also obtain an infinity of the two-integral 
odd DFs for this pair under the assumption of the rotation laws. It is worthwhile to mention that such equipotential extensions 
have been studied for the Jaffe (1983) spherical model (e.g. Jiang 2000; Jiang, Fang, Liu, Moss 2002; Jiang, Fang, Moss 2002) 
and for the spherical 7 model (Jiang & Ossipkov 2006). The 7 model was found independently by Kuzmin et al. (1986), Dehnen 
(1993), Saha (1993) and Tremaine et al. (1994). It is in fact an extension of the Jaffe model. The device of extending from 
the spherical models to the axisymmetric models was originally developed by Miyamoto and Nagai (1975) and independently 
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by Kutuzov and Ossipkov (1976). Other similar extensions can be also found in their further work (Nagai & Miyamoto 1976; 
Kutuzov 1989; Ossipkov 1997; Ossipkov & Jiang 2007a,b). 

Put V = — <&(R, z). Then, by (40), the mass density (41) can be expressed as a function p(ip, R) of tp and R, that is, it is 
of the following form 

p(V, R) = p(4>) + PW/X + M)/X 2 + pW/X 3 , (42) 
where and below, p(tp), p(ip), p(ip) and p(ip) are denoted by 
M rjw 4 + a 2 (3 + 2w)w 5 

^ ) = i^ ^TTp ' (43) 

_, M rjbw 3 (w + l) + 2a 2 b(3 + 2w)w 5 „, ,, M a 2 b 2 (3 + 2w)w 5 , M rja 2 bw 3 

pw = ^TIF —> pw = i^f (Ui)' ' pW = i^l^TTT (44) 

with u> = e rj ^^ GM ^ — 1. On the right side of (42), the first term is "spherical" and the other ones are axisymmetric like 
three types of the mass densities defined by (8), (9) and (10). The prolate Jaffe model can be hence regarded as a stellar 
system composed of a "spherical" component and three axisymmetric ones. Thus the DFs for the "spherical" and the three 
axisymmetric subsystems can be derived by using the formulae given in the previous sections and so the DF of the prolate 
Jaffe model can be expressed as a sum of the DFs of these subsystems. 

Throughout this paper, let H(x) denote the Heaviside step function, that is, H(x) — 1 as x > and H(x) = as x < 0. 
First, the two-integral even DF of the "spherical" subsystem can be expressed as 

v/Itt 2 de 2 J Q s/F=^' 
It can be also found from (29) that the two-integral even DF is of the form 



V8tt 2 dQ 2 Jo y/Q-ip 

corresponding to the second term on the right side of (42), where and below, Q = e — L 2 /(2a 2 ). 
We can first know from (44) that p(ip) can be expanded as a power series of the form 

p(i>) = J2 (47) 

71 — 5 

Then, by using (28), the even DF corresponding to the third term on the right side of (42) can be given by 



f+(e,L z ) 



(2^)3/2r(n-i/ 2 ) aFl I 1 ' 2- n '2'2^) as e>L z /(2a), 

(48) 



n— 5 



- + f {A) _1 (l, f ; " + % ) - e<Ll/(2a 2 ). 

n— 5 

Note that the power series (47) converges so fast that we can use a sum of N previous terms of (48) as a viable approximation 
of the even DF, where N is a natural number such that |ajv+5| is less than a small accuracy given in numerical calculation. 
Then this approximation, denoted by S+ (e, L z ), can be explicitly given by 

S N + {e,L z )=\ N+A ^ \ J (49) 

! - E (&) 2f\ (l,| ; n+i;2||i) as e < L 2 z /(2a 2 ). 

V. 71 = 5 

Hence the above series expression is computationally superior to the contour integral. 

By using (31), corresponding to the last term on the right side of (42), the two-integral even DF f + (Q,L z ) can be 
expressed as follows: 



f+(Q ' Lz) -7^ 



£ r pWdip l\ d 3 r p(ip)dip 

dQ 2 Jo a/Q - V ffl2 d Q 3 Jo y/Q-ip 



(50) 



Finally, by adding all the above four DFs (45), (46), (48) and (50), the two-integral even DF of the prolate Jaffe model 
with the mass density (42) can be obtained as a function /+(e, L z ) = f+(e) + f+(Q) + /+(e, L z ) + f+(Q, L z ). Thus, when 
e > L 2 /(2a 2 ), f+(s,L z ) is given by 

H(e) d 2 f e p(i,)dj, 1 d 2 f Q p(V>)# 
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+ 



g n r(n+l)e"- 3 / 2 / 3 _ 
2^ (2vr)3/2r(n- 1/2) 21 I ' 2 "'2' 2ea 2 / ^ V»Jr a 



Q p(ip)dt[> L 2 d 3 f Q p(ip)dt[> 



dQ 2 Jo VQ~i> a d Q 3 J o VQ-tP 



(51) 



when e < L 2 /(2a 2 ), /+(e, L z ) is of the form 



//(e) d 2 [ e ^ a n r(n + l)e"-^ / ^ V J. 3 1 2 £ a 2 \ 

f+{£ ' Lz) = 7^d^J 7^^2(2^/^ + 1/2) (.Wj 2Fl (, 1 '2' n+ 2'TTj- (52) 

Now let us consider the odd part of the DF for the prolate Jaffe model under the assumption of the rotation law. It can 
be known from (43) and (44) that the functions p(ip), p(ip) and p(tp) can be expanded as three power series: 

PW = *^2d n ip n , p(ip) = ^2 a n tp n , p(ip) = ^a„V n - (53) 

n— 4 n — 3 n— 3 

Similarly, by first inserting (47) and (53) into (42) and then using (39), we can obtain the two- integral odd DF 



a 2fc sgn(L z ) ^ . „_ 2 / L 2 Z \ a 2k 1 sgn(L,) ^ . n _ 2 ( 1 

/ - (E ' Lz) = 4^F(fc) L a " n!e H ( fc ' °> n - ^ W J + 4^ + 1/2) L a " n!£ H ( k > 2 



k,^,n- 1,1; 



2 ' ' ' 2ea 2 / 

sgn(L z ) ^ „ n-2-afi 1 i i L * \ . a 2fc ~ 3 sgn(L z ) ^ _ n _ 2 / 3 L 2 \ 

4.r(fcTTT^ a,in!£ M r 1 ' n " 1|1; 2^J + 4 ff r(fc + 3/2) 2. fl " nle e ( fe '2' n - 1 ' 1; w ) 



(54) 



for the prolate Jaffe model with the mass density (42), constrained by the rotation law of the form (vj>) = R 2k /(l + R 2 /a 2 ) k 
for fc > 0. In particular, when < [k] < 2, it can be deduced from (36) that the two-integral odd DF corresponding to the 
third term on the right side of (42) can be given by 

f (O T 1 = (k + l)sgn(L z )LfH(Q) rfM+ 3 [Q m ^ 
J " W ' zj 47r2% 2 r(fe + 2)r(l-fc )dQW+3i (Q - V) fc ° ' 
Hence, by (39) and (55), we can also obtain the odd DF 

- q 2fc sgn(L z ) ^ . n _ 2 / L 2 \ a 2fc - 1 sgn(L z ) ^ . „_ 2 ( \ L\ 

f - {£ > = 47 rr(fc) 2. EI I fc, 0, n — 1,1; tt^TJ ) + 4 ^(fc + 1/2) 2. a ^ e ( ' 2 ' H ' 2E0 2 

n=4 x 7 n=3 x 

(fc + l)sgn(/^f//(Q) #1+3 ^ WW + g^fej ,„ n!e "- 2 e (fc, n - 1, 1; (56) 



47r2 fc a 2 r(fc + 2)r(l - fco) dQW+3 7 (Q_^)fc 47rr(fc + 3/2) ^ ^'2' ' ' 2ea 2 / 

for the prolate Jaffe model with the mass density (42) under the above rotation law for [k] £ (0, 2] . 

All the above series are theoretically convergent to the corresponding Hunter-Qian contour integrals, that is, these series 
expressions of the DFs can be in fact used to recover Hunter and Qian's formulae of the DFs for the prolate models. The 
convergence of these series is fast enough that we can adopt a sum of finite previous terms of series as a viable approximation 
of the DFs for understanding of dynamical properties of the stellar system. Let 5o denote the desired small accuracy. Then 
this approximation of the DF (56), denoted by S[(e, Q,L Z ), can be given by 

a 2fc sgn(L z )^ x . __ 2 _ / n ^ , 1 L 2 Z \ , o*-V(i,)^ s . ~- 2mr A 1 , L 2 Z ' 



(57) 



S - (£ ' fti ^1^2« !£ H f fc, 0, n - 1, 1; — j + - ^ a n n! £ H ( fc, -, n ~ 1, 1; — 2 \ 

n=4: ^ ' n=3 * ' 

(fc + l)sgn(L z )L 2fc //(Q) dW+3 ^ a 2fc - 3 sgn(L z ) ^ 2 / 3 

+ 4^a 2 r(fc + 2)r(l-fc ) d Q W+ 3y o (Q_^)feo + 47rr(fc + 3/2) ^ "^'2'" ' ' 2ea 2 J 

where I is dependent of the convergence of the power series given by (53), for example, we can choose the natural number I to 
be such that \di+i\ + 5;+3| + a;+3 < So since we find that the convergence of these power series is fast enough. Fortunately, 
all the terms in these series can be expressed in terms of hypergeometric functions and it is easy to calculate numerically 
these hypergeometric functions. Therefore this approximation method is in general computationally preferable to the contour 
integral method, especially for the prolate Plummer model considered below. 

Of course, for the prolate Jaffe models, Jiang and Moss (2002) gave the complex contour integral expressions of the 
even DFs and obtained explicit formulae of both velocity dispersions and rotation curves. We can know from Hunter's (1977) 
formulae that these velocity dispersions for the prolate models can be expressed in terms of elementary functions of the two 
variables R and z. We can also find that the prolate models have anisotropic velocity distributions. More other properties of 
the prolate models have been given by Jiang and Moss (2002) including some figures such as the even DFs, velocity dispersions 
and rotational velocity. Their details are therefore not repeated here. 
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It is finally worth mentioning that a different way to recover Hunter and Qian's formulae of the DFs for the prolate 
models is to employ the Laplace-Mellin transformation given in Section 3.1. 



4.2 The prolate Plummer model 

The gravitational potential of the prolate Plummer model is given by 

*(R, z) = GM (58) 

y/WK 2 + a 2 + b) 2 + z 2 

where a, b, M and G are the same as in (40). As b goes to zero, (58) becomes a potential of the spherical Plummer (1911) 
model. The mass density corresponding to (58) can be calculated as 

M br 2 (X 2 + a 2 )+3a 2 X(X + b) 2 
P( R , z >~4^ ^3 ( 59 ) 

where X and r are the same as in (41). It can be found that the model defined by (58) and (59) is prolate and that for any 
given a > 0, the larger the value of 6, the more prolate this model. 

Similarly, by letting ip = — $(R, z) and using (58), the mass density (59) can be rewritten as a function p(ip, R) of ip and 
R, that is, it is of the following form 



6a 2 by 5 H 3 



4ttG 5 M 4 4nG 3 M 2 



1 3o 2 &V 1 b£ 1 

X 4ttG 5 M 4 X 2 4nG 3 M 2 X s ' 1 ' 



Take Q = e- L 2 /(2a 2 ). Then, by using (28), (29), (31) and (32), the two-integral even DF f + (e,L z ) corresponding to the 
mass density (60) can be obtained as follows. When Q < 0, 

, , r . ma 2 e 7 ' 2 H{e) QAaH 2 e 9 ' 2 / 3 11 2ea 2 \ 

/+(£ ' Lz) 7V8 W 3 G 5 M 4 ^ 2 iV8 7 r3 G 5 M 4 Li 2Fl ^' 2 ! Y ; TTj 5 (61) 
when Q > 0, 

96a 2 £ 7/2 Jj-(e) 2feQ 3/2 (96a 2 Q 2 + 7G 2 M 2 )H(Q) 
~ 7\/87r3G 5 Af 4 + 7\/87r3 a G 5 M 4 

96&V /2 / H L 2 \ bQ^ 2 (2a 2 Q-3L 2 z )H(Q) 

7^3 G 5 M 4 2 1 ^ ' 2'2'2ea2 y l + y87r3 ffl 5G3M2 ' 1 j 

Obviously, it is easier to calculate numerically the even DF given by (61) and (62) than the Hunter-Qian complex integral. 

Similarly, we can also give a simple expression of the odd DF for the prolate Plummer model constrained by the rotation 
curve. For example, if the rotation curve is assumed to be determined by (vj,) = R 2 /(l + R 2 /a 2 ), then, by (39), the odd part 
of the DF corresponding to the mass density (60) can be obtained as 



fJe,L z ) = a sgn(L z ) 



2* 2 &M± m °' 3 ' 1 '2rfj + ^/2G5M 4H I, 1 ' 2 ' 3 ' 11 W J 



3&e A 45a 2 6 2 e 3 L 2 \ be „„ / 3 



+ ,_ S /,^,,, ^ + n _ 9 ^ s ,^ H 1,1,3,1;^ + o_ S /,^,,, MI U.o.2,1 



4tt5/2G3M 2 V '2' ' ' 2ea 2 / 2tt 2 G 5 M 4 V ' ' ' ' 2ea 2 / 27r5/ 2 G 3 M 2 V '2' ' ' 2ea 



,2 



(63) 



All the functions H presented here can be expressed in terms of hypergeometric functions which can be calculated numerically 
with little additional effort. 



5 CONCLUSIONS 

It is recently shown that near the centre the potentials of the dark matter halos are approximately prolate on average 
(Hayashi, Navarro and Springel 2007). This implies that the study of prolate models is very significant for the investigation of 
the dark matter halos of the galaxies. To do this, it is a natural idea to explore different expressions of the two-integral DFs in 
axisymmetric stellar systems. These expressions presented here can be successfully applied into getting the even and odd parts 
of the two-integral DFs for both the prolate Jaffe and Plummer models. We can thus show that both the even and odd parts of 
the two-integral DFs for the prolate Jaffe model can be expressed as series convergent to the corresponding contour integrals 
obtained by use of Hunter and Qian's (1993) contour integral formulae. The convergence of these series is fast enough that we 
can use a sum of finite previous terms of series as a viable approximation of the two- integral DFs. All the terms of series can be 
expressed in terms of hypergeometric functions and it is easy to calculate numerically these hypergeometric functions. Hence 
this method is generally superior to the Hunter-Qian contour method, especially for the prolate Plummer model. Similar such 
expressions of series for the two-integral DFs have been considered for some stellar systems (e.g. Miyamoto 1971). 
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It is remarkable that the two-integral even DFs for some separable densities such as (9) and (10) can be expressed as a 
sum of finite real integrals which are easier to calculate numerically than their corresponding complex contour integrals. We 
can also find that the Hunter-Qian contour integral formula of the two-integral even DF for axisymmetric systems can be 
recovered by use of the Laplace- Mellin integral transformation originally developed by Dejonghe (1986). 

According to the Hunter (1977) formulae, a simple expression of the radial, vertical and rotational velocity dispersions 
can be directly given by using the axisymmetric potentials in the stellar systems. It can be also known that the radial velocity 
dispersion is equal to the vertical velocity dispersion for all the two-integral models in axisymmetric stellar systems. However, 
in real axisymmetric stellar systems, the velocity dispersion in the radial direction is not equal to the velocity dispersion in 
the vertical direction. This means that the DFs of the real systems must actually depend on three integrals of the motion. In 
axisymmetric systems, there are two classical isolating integrals, that is, the energy and the angular momentum with respect 
to the axis of symmetry. It can be also found by numerical calculation that the third isolating integral is respected by orbits in 
realistic axisymmetric galactic potentials. The third isolating integral cannot be generally expressed by using basic elementary 
functions of energy and angular momenta. But some three-integral models have been constructed for particular orbital families 
in flattened axisymmetric systems (Evans et al. 1997) and for separable axisymmetric Stackel potentials (Famaey et al. 2002). 

It is hence a very significant research field in both mathematics and astrophysics to consider the three-integral DFs for a 
good description of elliptical galaxies. Any three-integral model requires the existence of a global third isolating integral and 
its study is now still based on special potentials which allow the three-integral DFs to be obtained as tractable expressions. 
For example, Dejonghe and de Zeeuw (1988) first constructed the three-integral DFs for the Kuzmin-Kutuzov (1962) models 
since these models have Stackel's potentials. In general, either the prolate Jaffe potential or the prolate Plummer potential is 
not of the Stackel type that seems the most general of integrable potentials (Lynden-Bell 1962b) . The question addressed here 
is, is it possible to develop the three-integral DFs either for the prolate Jaffe model or for the prolate Plummer potential? This 
is an open problem now. To answer this question, we have first to see whether there exists a global third isolating integral 
for the prolate Jaffe model, which can be expressed in terms of energy and angular momenta. Once the third integral is 
found, a tractable expression of the three-integral DFs can be given as Dejonghe and de Zeeuw (1988) did. In principle, this 
is a reasonably clear procedure. However, so far it has only been carried through for special cases mentioned above; it seems 
difficult to obtain useful results for more general cases. 
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